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We investigate an attractive atomic Bose-Einstein condensate (BEC) trapped by a double-well 
potential in the axial direction and by a harmonic potential in the transverse directions. We obtain 
numerically, for the first time, a quantum phase diagram which includes all the three relevant phases 
of the system: Josephson, spontaneous symmetry breaking (SSB), and collapse. We consider also 
the coherent dynamics of the BEC and calculate the frequency of population-imbalance mode in 
the Josephson phase and in the SSB phase up to the collapse. We show that these phases can be 
observed by using ultracold vapors of 7 Li atoms in a magneto-optical trap. 
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In many experiments atomic Bose-Einstein conden- 
sates (BECs) are cigar-shaped due to a strong harmonic 
trapping potential in the cylindric radial plane; these 
BECs can be separated in two parts by means of a 
double- well potential in the cylindric axial direction 
This kind of geometry is the ideal setup to study the 
Josephson effect, a macroscopic coherent phenomenon 
which has been observed in systems as diverse as su- 
perconductors @, superfluid Helium 0] and, recently, 
also BECs in trapped ultracold atomic gases The 
observed coherent dynamics of the atomic BEC in the 
double-well potential (bosonic Josephson junction) [l|, 0| 
is efficiently described by nonlinear Josephson equations 
(JEs) 1, which are based on a two-mode approximation 
of the Gross-Pitaevskii equation (GPE) Q. These JEs 
are fully symmetric by changing the sign of the inter- 
atomic scattering length and do not predict the collapse 
of the BEC. The collapse of an attractive BEC of 7 Li 
atoms or 85 Rb atoms has been observed by various ex- 
perimental groups 0] and theoretically analyzed by many 
authors: in a single- wellpotential 0, in a potential with- 
out axial confinement in a toroidal confinement [Ioj |. 
in a double- well potential , and in a periodic potential 

In this paper, by correctly taking into account the di- 
mensional reduction of GPE from 3D to ID, i.e. by us- 
ing the so-called ID nonpolynomial Schrodinger equation 
(ID NPSE) [lH, we show that for an attractive BEC 
(negative inter-atomic scattering length) the JEs are not 
reliable in the presence of strong coupling. By numeri- 
cally solving the ID NPSE we obtain, for the first time, 
a quantum phase diagram of the three relevant regimes 
of the attractive BEC in a double-well: the Josephson 
phase, where the metastable state of lowest finite energy 
has a balanced population Q , the spontaneous symmetry 
breaking (SSB) phase, where the metastable state has an 
unbalanced population [l|, H[, and the phase of collapse, 



where the system reaches the collapsed ground-state with 
energy equal to minus infinity. Note that the problem of 
BEC collapse in an axial double well potential has been 
investigated by Sakellari, Proukakis, and Adams [ll|, but 
they have not derived the quantum phase diagram of the 
attractive BEC. Instead, very recently the collapse region 
in a quantum phase diagram has been obtained for a pai r 
of cigar-shaped traps coupled by tunneling of atoms [lj| . 
We also study the coherent dynamics of the system and 
calculate the frequency of the population imbalance both 
in the Josephson regime and in the SSB regime. In the 
SSB phase this frequency reaches its maximum value at 
the coupling strength where there is the collapse of the 
BEC. In addition, from the ID NPSE we obtain general- 
ized Josephson equations, which we call nonpolynomial 
Josephson equations (NPJEs) for the fractional imbal- 
ance and relative phase of the bosonic Josephson junc- 
tion. These new NPJEs reduce to the familiar JEs in the 
weak-coupling limit, but show a better agreement with 
the numerical results of the ID NPSE (and 3D GPE) for 
strong couplings (for both positive and negative scatter- 
ing length). Finally, we suggest that our predictions can 
be observed experimentally by using an ultracold vapor 
of 7 Li atoms and tuning the s-wave scattering length. 

Let us consider a dilute interacting BEC at zero tem- 
perature confined by a trapping potential Vt rap (r). This 
potential is taken to be the superposition of an isotropic 
harmonic confinement in the the transverse radial plane 
and a double- well potential Vow (%) in the axial direction 
x. Then, V tra p(r) is given by 



Vtrap(r) = V DW (x) + 



^P 2 



(1) 



where p is the cylindric radial coordinate, m is the mass of 
the atom, and wj_ is the trapping frequency in the radial 
plane. The macroscopic wave function ^(r,^ describing 
the above system with N atoms is governed by the 3D 
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where a s the s-wave boson-boson scattering length and 
*(r, t) is normalized to 1. The 3D GPE captures the 
main properties of collapse threshold and, as shown by 
using a reliable nonlocal potential, the collapsed state is 
actually a state of very high density which decays due to 
inelastic two- and three-body collisions [HI . By following 
Ref. jl3j ]. we choose the wave function ^(r,t) as the 
product of an axial complex wave function f(x,t) and 
a Gaussian transverse wave function of radial width a, 
where a depends on the axial wave function f(x,t), i.e. 
(7 = a(f(x,t)). By expressing lengths in units of a± = 
\/h/muj_\_, times in units of wj 1 , and energies in units of 
huj±, it is easy to show that the fields f(x,t) and a(x,t) 
satisfy the following equations [l3| 
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where T = 2a s N/a± and f(x,t) is normalized to 1. In- 
serting Eq. (01 into Eq. ([3]) one gets the so-called ID 
NPSE [Ljj], which is extremely accurate in reproducing 
the properties of the full 3D GPE with transverse har- 
monic confinement 1131. 
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FIG. 1: (color online). Axial probability density |/(a;)| 2 of 
the metastable attractive BEC in the symmetric double-well 
potential Vdw(x) where the two minima are at x = ±xo with 
xo — 2 and the energy height is Uo = 0.8. F = 2Na a /a± is the 
interaction strength. Results obtained by using ID NPSE, 
Eqs. ([3]) and Q. Length x in units of ai = y^h/muj±, 
density |/| 2 in units of aj 1 and energy in units of Tilo±. 

We have solved the ID NPSE by using a finite- 
difference Crank-Nicolson code with imaginary time [l6l . 
[l~7| to obtain the ground-state of BEC in the symmetric 
double-well trap. In the numerical analysis the double- 
well potential Vdw{%) is given by the combination of two 



Pbschl- Teller potentials with the energy barrier of height 
Uo = 0.8 huji_ and the local minima at —xq = —2<n_ and 
xo = 2a± (for details see Refs. [H, H^). It is important 
to stress that with T < the ground-state is always the 
collapsed state with energy equal to minus infinity. Thus, 
for r< Owe are actually looking for the metastable state 
of lowest finite energy. We find that this metastable state 
is symmetric for Tssb < T < (Joscphson phase), it has 
a broken symmetry for Tc < T < Tssb (SSB phase), 
and it becomes the collapsed ground-state for T < Tc 
(collapsed phase). In Fig. Q]we plot the axial probability 
density |/(a;)| 2 of the metastable state obtained by solv- 
ing the ID NPSE with imaginary time. We start with 
a slightly asymmetric initial condition and proceed up 
to the convergence to a stable configuration. The figure 
shows that for T = —0.2 the profile of the (meta-)stablc 
state is symmetric, while for T = —0.3 it is not. For 
r = —0.8 the BEC is practically localized only in the 
right well. In addition, we find that for T < —1.2 there 
is the collapse. 
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FIG. 2: (color online). Quantum phase diagram (Uo, \T\) of 
the attractive BEC in the symmetric double-well potential 
Vdw(x). The two minima are at x — ±xo with xo = 2. Uo is 
height the of the central barrier of Vdw(x) and V = 2Na s /a± 
is the interaction strength. There are three phases: Josephson 
(J), spontaneous symmetry breaking (SSB), and collapse (C). 
Solid lines are obtained with ID NPSE, Eqs. © and @. 
Filled circles and squares are obtained with 3D GPE, Eq. 
J2j. Energy Uo in units of hu}±. 



It is very interesting to analyze the quantum phases 
of the attractive BEC as a function of the height Uq of 
the energy barrier. We have performed a systematic in- 
vestigation by changing both Uq and T. The results are 
shown in Fig. [2] where we plot the quantum phase di- 
agram in the plane (Uq, \T\). To our knowledge, this is 
the first time that this kind of phase diagram is obtained 
for an attractive BEC in a symmetric double- well poten- 
tial. In the figure, the solid lines are obtained from the 
ID NPSE, Eqs. © and ©. The interaction strength 
|rsss| of the SSB transition strongly decreases by in- 
creasing the height Uq of the energy barrier, while the 
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critical strcgth \Tc\ to get the collapse slightly decreases 
by increasing Uq. To check the accuracy of the NPSE 
we have also solved the 3D GPE, Eq. by using a 

cylindric-symmetry finite-difference Crank-Nicolson code 
with imaginary time 17(. In Fig. [2] the collapse points 
predicted by 3D GPE are shown as filled circles, while the 
SSB points of 3D GPE arc filled squares. As expected 
[H|, the agreement between ID NPSE and 3D GPE is 
quite good. We stress, however, that both 3D GPE and 
ID NPSE are based on a zero-range inter-atomic poten- 
tial d, Il3j . A more accurate description of interaction 
might lead to a slightly different transition between the 
SSB and the collapse phase in Fig. [2] 

Let us now consider this two-mode approximation of 
the ID NPSE. Under the condition that the central bar- 
rier Uq of the double- well potential Vdw(x) is sufficiently 
high - corresponding to a weak link between its left ad 
right sides - the field f(x, t) can be decomposed by using 
the two-mode approximation 



f{x,t) = f L (t)<l> L (x) + f R (t)4> R ( X ) 



(5) 



The functions (/)l(x) and 4>r{x), which are orthonormal, 
are localized in the left and right well, respectively. We 
assume that the above functions are real and use the 
ansatz (|5|) in the NPSE. We multiply the resulting equa- 
tion by (j>ct(x) (a = L, R), and integrate over x. Then, by 
taking into account the overlaps between </> Q 's localized 
in the same well and neglecting those ones between <^ a 's 
localized in different wells, we obtain 
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where the parameters e, K and g are given by 
+°° r i d 2 



K = 



— oo 



dx(j) a (x) 
dx(frp(x) 



2dx 2 

iaP_ 

2dx 2 



V DW (x) 
Vdw{x) 



= r 



(7) 

Mx) (8) 
Mx) (9) 
(10) 



dx((f> a (x)) 4 . 

Finally, we write the time-dependent amplitudes f a {t) 
as f a (t) = y/N a (t) e te " m , with N a (t) being the fraction 
of bosons in the a-th well and 9 a (t) the corresponding 
phase. Then, Eqs. (|6|) and (J7J) give rise to the following 
system of coupled ordinary differential equations for the 
fractional imbalance z(t) = Ni,(t) — N R (t) (here N^t) + 
N R (t) = 1) and the relative phase 6(t) = 6 R (t) - 9 L (t): 
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VT^(4 + 3. 9 (1 -z)) 



(12) 



We call these equations " nonpolynomial Josephson equa- 
tions" (NPJEs) because they are derived from the NPSE. 
NP JEs describe the dynamics of the fractional imbalance 
z(t) and relative phase 8(t) of the bosonic Josephson 
junction taking into account transverse-size effects. It 
is clearly much easier to solve numerically these NPJEs 
than the full 3D GPE or the ID NPSE. When the cou- 
pling strength g is much smaller than one our NPJEs 
become 



i = -2 K V+ 

z 



e = 2K 



■. cos 8 + gz , 



(13) 
(14) 



which are the familiar Josephson equations (JEs) for a 
BEC found in Ref. Q. It is straightforward to verify 
that Eqs. (|13[) and (|14[) are invariant under the transfor- 
mations r — » — r and 8 — >• — 8 + it. Instead, Eqs. (fTl"j) 
and p2p do not exhibit this invariance. 

The stationary Josephson regime corresponds to the 
equilibrium points with 2 = and 8 = (balanced pop- 
ulation). Both JEs and NPJEs show that this station- 
ary Josephson phase exists only for Tssb < I\ where 
^ssb < 0. We have verified that the value of Tssb pre- 
dicted by JEs is always very close to the value obtained 
with NPJEs. The points below Tssb correspond to SSB 
phase. This phase, according to the JEs, exists for any 
Tssb < T; while, according to the NPJEs, the SSB phase 
does not exist anymore at the collapse strength Tc ■ Thus 
NPJEs predict a collapsed phase while JEs do not. 
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FIG. 3: (color online). Fractional imbalance z as a function of 
the time. Upper panel: V = -1.05. Middle panel: V = -0.25. 
Lower panel: T = 0.25. The solid lines are obtained with ID 
NPSE, Eqs. @ and ©. The dashed lines are obtained with 
NPJEs, Eqs. JTTJ and The dot-dashed line are obtained 

with JEs, Eqs. (fT3]) and (TTi)) . Initial conditions: 2(0) = 
0.2 and 8(0) — 0. Parameters of the double-well potential 
Vdw(x): energy barrier height Uo = 0.86 and location of the 
two minima at x = ±xq with xo = 2. Lengths in units of a±, 
time in units of ujJ 1 , energies in units of hu}±. 

To compare NPJEs with JEs we plot in Fig. [3] the 
population imbalance z(t) for three values of T, choosing 
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as initial conditions z(0) = 0.2 and 6(0) = 0. The figure 
shows that the NPJE curves (dashed lines) are always 
closer to the NPSE results (solid lines) than the JE ones 
(dot-dashed lines). Nevertheless, for a sufficiently strong 
(and negative) T the predictions NPJEs are no longer re- 
liable. Notice that while for T = 0.25 the system displays 
Josephson oscillations, i.e. coherent oscillations around 
z = 0, for T = 0.25 and T = -1.05 there are SSB oscilla- 
tions, i.e coherent oscillations around z = zssb ^ 0. 
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FIG. 4: (color online). Left panel: spontaneous symmetry 
breaking oscillation frequency ujssb around (z — zssb-,0 = 
0) vs. F. Right panel: Josephson oscillation frequency uij 
around (z = 0, 9 = 0) vs. T. Filled squares are obtained with 
ID NPSE, Eqs. (j3} and The dashed lines are obtained 
with NPJEs, Eqs. (JTTJ) and |12]). The dot-dashed line are 
obtained with JEs, Eqs. (|13[) an d dt - Parameters of the 
double well and units as in Fig. [3] 

We investigate in detail these coherent oscillations by 
looking for the stationary points and calculating the fre- 
quency of small oscillations around these points. In the 
case of NPJEs and JEs we diagonalize the Jacobian ma- 
trices associated, respectively, to the NPJEs (|13|) -(fT4 |) 
and JEs (JTT]) - (fT2")) . In general, the 2x2 Jacobian ma- 
trix has two complex eigenvalues, Ai^, and the stationary 
point is stable when Ai,2 = ^fiuj with ui > the frequency 
of stable oscillations. 

We consider first the Josephson regime and thus we 
study oscillations around the equilibrium points with 
3 = and 9 = (balanced population). Notice that 
for r = the oscillation frequency ujj reduces to the 
Rabi frequency, i.e. wj = 2K/fr. In the right panel of 
Fig. 0]we report the frequency wj of coherent oscillations 
around z = as a function of T. From the plots of Fig. 
2] one can see the differences between the behavior of u> j 
predicted by ID NPSE (dots), NPJEs (dashed line), and 
JEs (dot-dashed lines) . Among the three sets of data ob- 
viously NPSE ones are the more reliable. The softening 
of luj as T < approaches Tssb = —0.21 is reproduced 
extremely well by both NPJEs and JEs, while there are 



differences in corresponcncc of the hardening of OJj for 
large positive values of T. At T = 1.5 the relative error 
in the determination of uij between NPSE and NPJEs 
results is about 15%. 

We now analyze the oscillations in the SSB regime. As 
previously stressed to this regime is associated a symme- 
try breaking of the fractional population imbalance z, i.e 
the stationary configuration has 8 = but z = z$sb =h 0. 
In the left panel of Fig. [4j we plot the SSB oscillation 
frequency cos SB around the stationary z$sb ^ (with 
6(0) = 0) as a function of T. From these plots one can 
see that when T is small enough the behavior of ujssb pre- 
dicted by NPJEs and JEs are quite similar, but NPJEs 
results arc slightly better than JEs ones, cos SB is equal 
to zero at Tssb and it increases by decreasing T < 0. 
As r approaches the collapse strength Tc = —1.1 of ID 
NPSE the relative error in the determination of ujssb be- 
tween NPSE and NPJEs results becomes quite large: it 
is about 50% at V = —1.1. We have verified that, in con- 
trast with JEs, NPJEs predict the BEC collapse but at a 
critical strength much smaller (in modulus much larger) 
than the one obtained by using the NPSE. 

It is important to stress that it is possible to achieve 
the so-called "self-trapping regime" in correspondence of 
initial conditions (z(0),6(0)) which are not stationary 
points of JEs [f|. This regime is characterized by popu- 
lation imbalance (z(t) ^ 0) and a running phase during 
the time evolution. By solving our NPJEs we find this 
dynamical self-trapping for r c < T < FrgT,-) < or 
T > T(st,+) > with the thresholds F(5 T =F ) depend- 
ing on the initial conditions z(0) and 6(0). Note that by 
solving JEs one finds r c = and also Tst,- = — Pst,+- 

In conclusion, we observe that the results obtained so 
far, can be used to describe concrete systems. For in- 
stance, by considering an attractive Bose-Einstein con- 
densate made of 7 Li atoms, and choosing the transverse 
confining frequency as co±_ ~ 2-7T x 100 Hz, we have 
a typical value of the transverse length a± ~ 4 /xm, 
while the parameters of the double-well potential read: 
U ^ 6 • 10~ 32 J and x ^ 8 fim [1 [20[. The natural 
scattering length of 7 Li atoms is a s = —1.45 nm but 
it can be modified with an external constant magnetic 
field by means of a Feshbach resonance 0- Working 
with N ~ 10 4 condensed atoms in the trap, it is possible 
to observe experimentally the behavior of Josephson fre- 
quency ujj and of the SSB frequency ujssb by tuning the 
scattering length a s from positive values to the collapse 
point at a s = T c a±/(2N) ~ —0.2 nm. Finally, we stress 
that with the above values the condition \g\/K N 2 is 
fully satisfied, and the system is always in the coherent 
regime [IE [El. 

The authors thank Flavio Toigo for enlightening dis- 
cussions. 
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